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ABSTRACT 


The  effects  of  trim  on  stability  of  motion  during  depth 
control  of  submersible  vehicles  are  analysed.  Full  state 
feedback  control  is  used  to  provide  stable  response  in  the  dive 
plane,  and  feedforward  control  is  used  to  ensure  steady  state 
accuracy.  A  comple  Sr '  of  stability  maps  is  generated  for 
various  values  of  meta  »rtric  height,  longititunal  center  of 
gravity/center  of  buoyancy  separation,  forward  speed,  and 
control  law  time  constant.  The  results  clearly  indicate  ranges 
of  parameters  that  should  be  chosen  in  design  and  operation  of 
a  given  vehicle. 
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I .  INTRODUCTION 

The  fundamental  goal  of  submarine  control  is  to  reach  and 
maintain  ordered  depth.  Experimental  designs  involve  expensive 
model  testing  such  as  Oarpa  Suboff  Model  (DTRC  Model  5470) 
[Ref.  6].  Much  research  has  been  done  in  depth  control  of 
submarines  [Ref.  3,5].  Our  goal  is  to  develop  an  analytic 
method  to  determine  the  stability  properties  of  a  design. 

The  stability  of  a  design  will  have  a  significant  impact 
on  its  responsiveness.  A  vehicle  with  a  large  margin  to 
instability  will  not  be  very  responsive.  The  problem  becomes 
one  of  determining  how  close  to  the  margins  we  can  get  without 
a  total  loss  of  stability.  By  expanding  the  scope  of  our 
research  to  include  nonlinear  terms  we  are  able  to  define  the 
limits  of  stability  and  therefore  margins. 

Previous  studies  analyzed  stability  properties  of  the 
system,  specially  static  bifurcations  [Ref.  2]  and 
bifurcations  to  periodic  solutions  [Ref.  1].  The  latter  study 
which  is  used  as  a  basis  for  this  work,  was  restricted  to 
level,  zero  trim  flight  paths. 

The  purpose  of  this  thesis  is  to  develop  a  program  for 
finding  the  limits  of  stability  for  an  out  of  trim  submarine 
at  moderate  and  high  speeds.  These  limits  are  mapped  using  a 
Hopf  bifurcation  analysis  program  included  in  the  Appendix. 
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II.  EQUATIONS  OF  NOTION 

The  motion  of  the  submersible  in  the  vertical  plane  can  be 
modeled  by  four  coupled  nonlinear  differential  equations  for 
pitch  rate  (q),  heave  velocity  (w) ,  pitch  angle  (0)  and  heave 
(z).  With  a  body  fixed  coordinate  frame  at  the  vehicle's 
geometric  center  ,  we  can  express  Newton's  equations  of  motion 
as 


mfw-Uq  -Zffl2-  X(fl)  =  Zfi  +  Zji>  +Z6&b  +  Zfit  ♦  Z„w  +  Zp 


(w-xq)3  ^ 


(2.1) 


Iy-mxa(*  -Uq)  -  zawqm  -  -x^Bcos  0  -  zmBztn  0 
Mi*b  +  +  +  +  “P/  CDb <*>  *** 


(2.2) 


6  =  q 


(2.3) 


i  =  -  Uzin  0  +  w  cos  0 


(2.4) 
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Equations  2.1  through  2.2  can  be  written  in  a  more  compact 


form  as, 

i*  =  ajjUw*  aJ2Uq  +  a^z^tin  0  +  a^x^coa  0  ♦  bjU2bt 

b2U26b+  dw(w,q)  *  Cjfw.q) 


(2.5) 


q  =  a2JUw  +  a^Uq*  a^g^tbt  0  +  a2Jxagcoi  0  +  bjU2 4f 
b2U26b+  dq(w,q)  +  c3(»,q) 


(2.6) 


where , 

Dv  =  (m  -ZJQy-MJ-  (mxQ+  ZQ)(mxa  +  MJ 

°uD v  =  (Z„-2CDE0Utan  e^dy-Ai^)  ♦  (mxQ+ZJ(Mw+  IC^^Um  0^ 

anD*  =  (m+zq+  2 CdEj Utem  0 o)(Iy-M^  + 

(mxQ+  Zq)(Mq-  mx q  - mz QU tan  0Q-  2CDE2Uum  0^ 

=  2CDEjUtan  O^fm  -  Z^)*  (mxq  +  Mty(Zw-  2CDE0Utan  0^ 

"22-^,  “  (i/€-«*o-ii«*oVtm0o- 2C^2C/tei0o) 

(mxg  +  A/^Hm  +  Z?+  2  tan  0O) 

(2.7) 

m  -(m*  Z*)B 

a13Dv  =  *(»«o+ VB 

(2y~Ai^Zb~  (q+  Zj)Mb 
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b2Dv  =  (m  -  ZJMt+  (mx0*MJZt 


dq(w,q)Dv  =  (m  -  ZJlq+  (mxq  +  M+)I„ 

dw(w,q)Dv  =  (Iy-M4)Iw-  (mxa+  Z^Iq 

ct(w,q)D w  =  (Iy-M^mzQq2-  (mx0  +  Z^mzQwq 

c2(w.q)Dv  =  -  (m  -  Z+)mzQwq-  (mx  Q  *  M+)  m  zQq  2 

In  these  equations  the  submersible  is  assumed  to  be 
neutrally  buoyant  (W«B) ,  and  statically  stable  (zG  >  zB) .  Here 
we  can  assume  zB  to  be  zero,  hence  ZGB  ~  ZG. 

At  steady  state  the  cross  flow  drag  integral  terms  IH  and 
I,  have  the  form, 

/*  =  -CDw\w\fb(k)dc  Ip  =  CDw\w\J  b(x)xdx  (2.8) 

From  equation  (2.3)  it  is  seen  that  w  is  equal  to  tan90  at 
steady  state.  The  Jb(x)dx  term  is  computed  numerically  for  the 
SUBOFF  model  as  E0,  and  /b(x)xdx  term  as  Eu  Therefore,  the 
cross  flow  drag  terms  become, 

I„  =*  -CB  w  /W-/  E0  IP  -  CDw  /W  Ej  (2.9) 
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Because  we  have  two  rudders  at  the  bow  and  the  stern,  our 
system  of  equations  is  multi-input.  To  reduce  this  system  into 
a  single  input  system  the  linear  combination  of  the  control 
inputs  will  be  modified  into  the  following  form, 

5  -  5.  ,  5b  «  ad.  (2.10) 
where  a  is  the  control  surface  coordination  gain.  The  value 
of  a.  ranges  from  -1  to  1.  The  selection  of  the  value  of  a  will 
allow  the  planes  to  operate  as  desired  for  the  particular 
maneuvering  condition,  i.e.,  a  **  0  for  no  bow  plane  control, 
a  =  -i  for  bow  plane  and  r*e rn  plane  opposed  to  each  other, 
yielding  the  maximum  pitch  moment,  and  a  »  1  for  bow  and  stern 
plane  control  in  the  same  direction,  yielding  the  maximum 
heave  force. 
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FIGURE  2.1  Geometric  representations  of  the  basic 
definitions. 
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Ill 


CONTROL  LAX 


A.  INTRODUCTION 


The  control  design  problem  can  be  expressed  in  state  space 
as  follows, 


i  =  Ax  +  Bb 


(3.1) 


where  the  state  vector  is 


x  = 


(3.2) 


Equation  3.1  in  our  case  is. 


6 

0 

0  l  0 

0 

* 

a13sGB  ~  V  2*l 

ajju-hju’k,  ajfi-tjuh,  -bju2k4 

w 

4 

a23ZGB~b? 

a31u-bji3k2  arf-bjthj  -bjt3k4 

t 

-u 

1  0  0 

Z' 

(3.3) 


Our  aim  is  to  find  a  controller  which  will  assure  us  a 
stable  closed  loop  system.  The  only  control  input  is  the  dive 
plane  angle ,  5 . 
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B.  FEEDBACK  CONTROL 
1 .  Pol*  Placement 

The  full  state  feedback  controller  is  a  linear  function  of 
the  states  and  has  the  form, 

6  -  -  K  x  (3.4) 

where  K  is  the  vector  of  feedback  gains  which  are  to  be 
determined  in  order  to  give  the  desired  closed  loop  system 
dynamics.  Substituting  equation  2.12  into  2.10  yields, 

i  =  (A  -  BK  )  x  (3.5) 

The  feedback  gains  K  must  be  chosen  such  that  A  -  BK  has  the 
desired  eigenvalues.  The  actual  characteristic  equation  of  the 
closed  loop  system  is  given  by, 

det(A  -  BK  -  al)  -  0  (3.6) 

The  required  values  of  K  are  obtained  by  matching  coefficients 
in  the  two  polynomials  of  the  actual  and  the  desired 
characteristic  equations.  Equation  3.5  becomes. 


6 

o  oio 

0 

0 

w 

alfOB  alP  alf  0 

w 

+ 

bjU  3 

4 

*23*08  a2JU  *2?  0 

f 

b2u3 

t 

-II  1  0  0 

z 

0 

(3.7) 


The  characteristic  equation  of  the  closed  loop  system  is. 


-s  0  10 

al?OB-bp2kl  alJu-bJu3k2't  aJ^t~bJu2k3  ~bp2k4 

a23taB~b^2kI  a2tU~b^2k2  b#3k4 

-II  1  0  -t 


(3.8) 


which  reduces  to. 


4  *  *  (J^2k2  *  ^ 3k3  ~  & j) S  ^  +  (-Bjkj  -  &2k2  ~  ^3k3~  ^4k4  "  + 

(-Clkl-C2k2-C3k3-Es)t*(-D1k4-D2k4)  =  0 


(3.9) 


where , 


A2 

* 

-B4 

»  b 2 

u2 

a3 

- 

*  b2 

u2 

B, 

(*22 

bi  - 

*i2  b2)  u3 

b3 

* 

C \  - 

(*n 

b2  -  aai  b2) 

u3 

c. 

3 

D1  = 

(  *23 

*>i  ~*13  b2  ) 

ZGB 

u2 

c4 

3 

(*22 

bl  + 

ba  -  a12  b2) 

u3 

d2 

« 

(*11 

b2  - 

*21  b2)  u 4 

(*11 

+  a22 

)  u 

e2 

3 

*23ZGB  +  (*12*21  ~  *11  & 

*22) 

u2 

b3 

" 

(  *13 

*21  ~  *11  * 23 )  ZGB 

u 

Now, 

let' 

s  assume  that  we 

want 

(3.10) 


at  -pj,  -pa,  -p,,  -p4  to  have  the  desired  system  response.  Then 
the  desired  characteristic  equation  is, 

432  h  (3*11) 

t*  +  ajt  *  a3*  +  *3*  *  *4  *  0  '  • 


9 


where , 


- 

Pi  +  Pi  + 

Pa  +  P* 

a2  - 

Pi  Pa  +  Pi 

Pi  +Pi  Pi  +  Pi  Pi  + 

Pi  Pi  +  Pi  Pi 

(3 

.12) 

a3  - 

Pi  Pz  Pi  + 

Pi  Pa  Pi  +  Pi  Pi  Pi 

+  Pa  Pa  Pi 

a4  - 

Pi  Pa  Pa  Pi 

The 

feedback 

gains  can  now 

be  computed  by  equ 

g 

the 

coefficients  of  equation  3.9  and  3.11, 

A2  k2  ^  Aj  iCj  *  ™  Oj  • 

Bj  k2  +  B2  k2  +  k3  +  Bt  k4  -  a2  +  E2  (3.13) 

Cj  kt  +  C2  k3  +  C<  k4  •  a3  +  Ej 
(D2  +  D2)  k4  - 

We  established  the  method  for  placing  the  poles  of  the 
system,  but  we  also  need  to  know  the  desired  locations  of  the 
poles . 

2.  Pole  Location  Selection 

In  a  typical  second  order  system  control  law  design  , 
transient  response  specifications  are  given.  This  results  in 
an  allowable  region  in  the  s-plane  where  the  desired  location 
of  the  poles  can  be  obtained.  For  higher  order  systems  the 
concept  of  dominant  roots  can  be  employed.  In  selecting  poles 
the  physics  of  the  system  must  be  considered.  If  the  poles  are 
too  negative,  a  small  time  constant  will  result,  and  the 
system  may  not  be  able  to  react  that  fast.  If  we  use  big  gains 
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K,  this  also  means  that  the  control  effort  will  be  large.  In 
practice  there  are  limits  based  on  actuator  size  and 
saturation. 

Considering  the  control  law  design  to  stabilize  the 
submarine  to  a  level  flight  path  at  0  ■  0  it  is  required  that 
the  submarine  return  to  the  level  flight,  after  some  small 
disturbances  in  6  or  z,  within  the  time  it  takes  for  the 
vehicle  to  travel  three  ship  lengths.  Since  our  model  is  14 
feet  long  and  its  velocity  is  5  feet/sec.,  the  required 
recovery  time  is  about  10  seconds.  This  means  the  time 
constant  is  3  seconds  and  the  closed  loop  poles  should  be 
placed  at  approximately  -0.3. 

After  placing  the  poles  using  equation  3.13  the 
control  law  is  found  to  be, 

5  -  0.9917  9  +  0.8333  w  +  0.6026  q  -  0.0351  z  (3.14) 

C.  FEEDFORWARD  COHTROL 

The  previous  discussion  on  feedback  controller  assures 
closed  loop  stability,  but  it  acts  as  a  regulator  in  other 
words  takes  all  the  states  to  zero  values.  If  we  have  constant 
disturbances  or  we  want  to  track  some  reference  value  other 
than  zero  we  can  not  do  this  with  state  feedback  alone. 
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FIGURE  3.1  Feedback  and  feedforward  control  application  to 
our  linearized  model 


In  the  case  of  non- zero  set  points  or  constant 
disturbances  we  again  need  to  have  the  exact  same  full  state 
feedback  to  ensure  closed  loop  stability.  But  we  also  need  to 
introduce  an  additional  term  to  our  controller  in  the  form 

u  -  -  K  x  +  k0  (3.15) 
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where  kg  is  the  constant  feedforward  tern,  kg  is  given  as 
(Ref.  4 ] , 

k0  -  Hc'l(0)  x0  (3.16) 

where  Xg  is  the  reference  values  of  the  states  and  Hg'^s)  is 
the  closed  loop  transfer  function, 

He'l(a)  -  C  (al  -  A  +  BK)'1  B  (3.17) 

Another  way  of  getting  kg  is  looking  at  the  steady  state 
equations  of  notion.  In  steady  state  all  the  tine  derivatives 
in  equations  2*1  through  2.4  go  to  zero  and  we  have. 


w  =  tart  0fl 
Z|i  +  Zj*  -  Ig  =  0 


(3.18) 

(3.19) 


-(*0~  xB>Bcog  Q0~  *<aBainB0  +  *  Uww  *  !P  =  0 


(3.20) 


If  the  equation  3.19  is  nultiplied  with  Mg  and  equation  3.20 
with  Z6  and  set  equal  to  each  other  and  plug  in  the  equation 
3.18,  we  have  an  equation  depending  only  on  6. 

<ZJih-MJh)kmS0+xQBZhco,*0+z0BZtMtn*0+ 

CD(E^b  -  EjZt)  tan  0,  \Um  O0|  1 4  ’ '  ' 

Where  E0  and  E1  are  the  integral  terns  computed  numerically  for 
the  Suboff  model.  The  steady  state  value  of  0O  is 
found  from  equation  3.20  by  using  a  Newton-Raphson  method 
in  Bifurl  program  in  the  Appendix.  Then  we  can  get  5  from 
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equation  3.19, 


g  _  *H~  ^ww  _  ~CDE9tm  dg  \ltm  69l  -  Zww  (3.22) 

z.  z* 

After  getting  5,  we  easily  find  k,,  from  equation  3.15, 

k0  -  5  +  JCx  (3.23) 

A  plot  of  the  steady  state  angle  60  versus  x,  for  different 
values  of  zQ  is  shown  in  Figure  3.2. 
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THETA  (degrees) 


XBG  (%  of  L) 


FIGURE  3.2  Steady  state  pitch  angle  Q0  as  metacentric 
height  varies 
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XV.  BIFURCATION  ANALYSIS 


A.  STABILITY 

The  nonlinear  equations  of  motion  in  the  dive  plane  2 . 1 
through  2.4  can  be  expressed  in  a  compact  form  as  follows, 

j t=f(x) 

Where  x  is  the  state  variable  vector  x-[0,w,q,z].  We  know 
that  the  equilibrium  points,  x0  of  the  system  are  defined  by. 


f(x0)  -  0  (4.2) 

This  is  a  nonlinear  system  of  algebraic  equations  and  it 
may  have  multiple  solutions  in  x„,  which  means  that  the 
nonlinear  system  may  have  more  than  one  positions  of  static 
equilibrium.  If  we  pick  one  equilibrium,  x0  we  can  establish 
its  stability  properties  by  linearization.  The  linearized 
system  becomes. 


t  *Ax 


(4.3) 


where  A  is  the  Jacobian  matrix  of  f  (x)  evaluated  at  Xo, 


dx* 


(4.4) 
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and  the  state  has  been  defined  to  designate  small  deviations 
from  the  equilibrium  x0, 

x  x-x0  (4.5) 

In  system  dynamics  as  long  as  all  eigenvalues  of  A  have 
negative  real  parts,  we  know  that  the  linear  system  will  be 
stable.  This  means  that  the  equilibrium  x0  will  be  stable  for 
the  nonlinear  system  as  well.  This  is  in  fact  Lyapunov's 
linearization  technique. 

B.  BIFURCATION 

Values  of  the  nonlinear  system  parameters  at  which  the 
qualitative  nature  of  the  system's  motion  changes  are  known  as 
critical  or  bifurcation  values.  The  phenomena  of  bifurcation, 
i.e.,  quantitative  change  of  parameters  leading  to  qualitative 
change  of  system  properties,  is  the  topic  of  bifurcation 
theory.  Euler  buckling  (Pitchfork  bifurcation),  limit  cycles 
(Hopf  bifurcation)  are  common  examples  of  bifurcation. 

Classical  definition  of  stability  states,  that  the  real 
part  of  all  the  eigenvalues  of  the  system  must  be  negative. 
Therefore,  our  initial  investigations  into  the  stability  of 
the  SUBOFF  model  was  to  find  those  eigenvalues  whose  real 
parts  cross  the  imaginary  axis.  We  used  the  bifurcation 
analysis  program,  included  in  the  Appendix,  to  calculate  the 
eigenvalues  of  the  system. 
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By  linearizing  the  equations  of  motion,  equations  2.1 
through  2.4  ,  the  state  space  equations  of  the  dynamic  system 
can  be  written  in  the  form, 

i=Ax+Bu  (4*5) 


where , 


u  -  -  K  x  (4.6) 

and  K  is  the  vector  of  controller  gains,  as  calculated  by  pole 
placement  in  equations  3.13.  The  eigenvalues  of  the  system  are 
found  by  solving, 

det|A-BK-sI|  -0  (4.7) 

In  the  bifurcation  program  a  pseudo-root  locus  method  is 
employed  where  the  time  constant,  Tc,  is  fixed.  The  constant 
Tc  fixes  to  placement  of  the  system  poles  at  a  given  nominal 
forward  speed  U0  and  then  the  model  speed,  U,  is  varied 
incrementally  with  the  system  eigenvalues  calculated  at  each 
speed  increment.  When  the  real  part  of  an  eigenvalue  changes 
sign  between  the  limits  of  a  speed  increment  a  bisection 
method  is  employed  to  find  the  speed  where  the  real  part  of 
the  eigenvalue  equals  zero. 

For  each  point  where  the  real  part  of  an  eigenvalue 
crosses  the  imaginary  axis  the  associated  Tc  and  U  are  plotted 
on  a  bifurcation  map.  This  map  delineates  the  regions  of 
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classical  stability  (all  eigenvalues  on  the  left  hand  plane) 
from  the  regions  of  instability.  A  family  of  bifurcation  maps 
were  generated  by  varying  nominal  speed,  U0,  initial 
stability,  z0B,  and  longitudinal  center  of  gravity /buoyancy 
separation,  xOT  of  the  submersible. 

Figure  (4.1)  shows  a  typical  bifurcation  map  with  its  five 
distinct  regions  [Ref.  1].  Region  I  is  the  area  of  classical 
stability.  In  region  II  there  is  one  real  positive  eigenvalue 
which  is  indicative  of  a  pitchfork  bifurcation.  Pitchfork 
bifurcations  of  this  model  were  previously  examined  by  Reidel 
[Ref.  1],  Regions  III, IV,  and  V  have  at  least  one  pair  of 
complex  conjugate  eigenvalues  with  a  positive  real  part.  This 
would  indicate  that  there  should  be  an  unstable  oscillatory 
behavior  for  the  model. 

C.  RESULTS  AND  DISCUSSION 

The  classical  stability  region  in  bifurcation  maps  lies 
between  pitchfork  and  Hopf  bifurcation  boundaries.  The  limits 
or  parameters  must  be  defined  for  the  system  designer  prior  to 
starting  the  design.  By  maximizing  the  region  of  stability  we 
can  give  the  decigner  the  most  leeway  in  his  work.  There  are 
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FIGURE  4.1  A  typical  bifurcation  map  showing  the  five 
distinct  regions 
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two  parameters  that  we  used  to  change  the  bifurcation  maps, 
the  longitudinal  separation  between  center  of  gravity  and 
center  of  buoyancy,  x0B  and  the  initial  stability,  z^. 

First  we  look  at  the  change  in  x^.  In  figures  4.2  through 
4.7  we  plotted  bifurcation  curves  for  different  initial 
stabilities  as  x<,B  varies.  We  can  see  that  as  x0B  increases 
the  Hopf  bifurcation  branches  HI  an  H2  move  towards  higher 
speeds  and  time  constants  and  thus  increasing  the  stability 
area.  The  H3  branch  however  remains  constant. 

The  other  important  point  that  we  observed  is  that  the 
system  becomes  unstable  at  nominal  speed  at  higher  time 
constants.  This  is  unexpected  because  we  are  designing  around 
our  nominal  speed.  A  more  careful  examination  in  the  trimmed 
case  shows  that  the  actual  forward  velocity  becomes,  /u2+w2. 
Therefore  the  system  may  become  stable  at  a  value  of  u  other 
than  nominal. 

The  next  parameter  we  examined  was  the  initial  stability, 
zGB.  Figures  4.8  through  4.14  show  the  effect  of  varying  za 
from  .2  to  .4  ft  for  different  x^  values.  The  H3  branch 
remains  constant  while  the  upper  speed  H2  branch  moves  down 
effectively  decreasing  the  area  of  stability.  The  low  speed 
curve  HI  moves  upwards  and  increases  the  low  speed  area  of 
stability. 
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VELOCITY  (NONDIMENSIONAL) 


FIGURE  4.3  Bifurcation  map  as  xg  changes  between  xg*0  and 
xg« . 2 ,  zg» . 2 . 
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Bifurcation  map  as  xg  changes  between  xg»0  and 
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FIGURE  4.7  Bifurcation  map  as  xg  changes  between  xg»0  and 
xg« . 3 ,  zgm . 4 . 
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FORWARD  VELOCITY  (NONDIMENSIONAL) 


FIGURE  4.8  The  effects  of  changing  zg  on  the  bifurcation 


maps,  xg=* . 3 
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The  effects  of  changing  zg  on  the  bifurcation 
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FIGURE  4.10  The  effects  of  changing  zg  on  the  bifurcation 
maps,  xg*.l. 
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TIME  CONSTANT 


FIGURE  4.11  The  effects  of  changing  zg  on  the  bifurcation 
maps,  xg*0. 
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FIGURE  4.12  The  effects  of  changing  zg  on  the  bifurcation 
maps ,  xg*- . 1 . 
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FIGURE  4.13  The  effects  of  changing  zg  on  the  bifurcation 
maps ,  xg»- . 2 . 
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FORWARD  VELOCITY  (NONDIMENSIONAL) 


FIGURE  4.14  The  effects  of  changing  zg  on  the  bifurcation 
mapa ,  xg=- . 3 . 
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V.  CONCLUSIONS  AMD  RECOMMENDATIONS 


A.  CONCLUSIONS 

Hopf  bifurcation  analysis  is  a  very  useful  design  tool  in 
the  design  and  evaluation  phase.  Hopf  bifurcation  analysis  and 
an  identification  program  that  can  evaluate  the  hydrodynamic 
coefficients  for  the  submersible  vehicle  will  be  very  useful 
and  save  money  and  time  by  reducing  the  amount  of  model 
testing.  An  effective  set  of  control  system  parameters  can  be 
generated  in  this  process  that  will  be  optimal  for  the  final 
design  of  the  submersible. 

This  type  of  analysis  can  set  the  limits  of  the  ranges  of 
important  parameters  such  as  metacentric  height  and 
longitudinal  separation  of  buoyancy /gravity  centers.  As  we 
have  seen  changes  in  these  two  parameters  can  have  dramatic 
effects  on  stability.  It  was  found  that  the  moderate  speed 
region  of  stability  increases  with  increasing  metacentric 
height.  The  same  is  not  true,  however,  for  high  speeds.  The 
longitudinal  separation  of  center  of  gravity/buoyancy  can  have 
a  profound  effect  on  stability.  It  was  found  that  the  vehicle 
may  be  unstable  even  at  nominal  speed.  This  was  attributed  to 
the  fact  that  at  high  trim  angles,  the  feedback  gains  which 
are  computed  at  zero  trim,  can  no  longer  guarantee  stability. 
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B.  RECOMMENDATIONS 


The  bifurcation  analysis  program  should  be  expanded  to 
evaluate  the  performance  of  the  submarine  including  effects  of 
external  forces  such  as  wave  effects,  currents,  and  free 
surface  effects. 
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APPENDIX  -  BIFURCATION  ANALYSIS  PROGRAM 


C  PROGRAM  BIFUR1  FOR 
C  BIFURCATION  ANALYSIS 
C  PARAMETERS  ARE:  TC  VS.  U 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  K1  JC2JO^C4JL,MQDOT,MWDOT,MQ,MW,MDS,MDB^!D, 

&  MASS.I  Y,P  1  ,P2,P3,P4,XGB,ZGB 

DIMENSION  A(4,4),FV1(4),IV1(4)7ZZ(4,4),WR(4),WI(4),XL(25), 

&  BR(25),  VEC0(25),  VEC 1  (25),  VEC2(25) 

COMMON  P1,P2,P3,P4 
C 

OPEN  (1 1  ,nLE- BIF 1  .RES,,STATUS=rNEW) 

OPEN  ( 1 2,FELE- BIF2.  RES',  ST  ATU  S='NEW’) 

OPEN  ( 1 3 ,FILE-BIF3  RES’, ST ATUS='NEW) 

C  NUMERIC  INFO  OF  DARPA  SUBOFF  MODEL 
WEIGHT=1 556.2363 
BUO  =1556.2363 
L  =13.9792 
IY  =561.32 
G  =32.2 

MASS  =WEIGHT/G 
RHO  =1.94 
XB  =0.0 
ZB  =0.0 
CD  =0.4 


37 


CD  =0.5*CD*RHO 


C 

WRITE  (*,*)  'ENTER  MIN,  MAX,  AND  INCREMENTS  IN  Tc  (nondim)' 
READ  (*,*)  TCMIN.TCMAXJTC 

WRITE  (*,*)  'ENTER  MIN,  MAX,  AND  INCREMENTS  IN  U  (nondim)’ 
READ  (*,*)  UMIN,UMAX,IU 
C  WRITE  (*,*)  'ENTER  NOMINAL  SPEED' 

C  READ  (V)UO 

WRITE  (♦,♦)  'ENTER  XG  AND  ZG' 

READ  (V)  XG,ZG 
U0=9 
C 

ZGB=ZG-ZB 
XGB=XG-XB 
TCMIN=TCMIN*L/UO 
TCMAX=TCMAX*L/UO 
UMIN  =UMIN*U0 
UMAX  =UMAX*U0 
C  HYDRODYNAMIC  COEFFICIENTS 
ZQDOT=-6.3300E-04*0.5*RHO*L**4 
ZWDOT=-1.4529E-02*0.5*RHO*L**3 
ZQ  =  7.5450E-03*0.5*RHO*L**3 
ZW  =-1.3910E-02*0.5*RHO*L**2 
ZDS  =-5 .6030E-03  *0.5  *RHO*L*  *2 
ZDB  =-5 .6030E-03  *0.5  *RHO*L*  *2 
MQDOT=-8.8000E-04*0.5*RHO*L**5 
MWDOT=-5.6100E-04*0.5*RHO*L**4 
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MQ  =-3.7020E-03*0.5*RHO*L**4 
MW  =  1.0324E-02*0.5*RHO*L**3 
MDS  =S-2.4090E-03*0.5*RHO*L**3 
MDB  =  2.4090E-03*0.5*RHO*L"3 


XL(1>=0.0 

XL(2)=0. 1 

XL(3)=0.2 

XL(4)=0.3 

XL(5)=0.4 

XL(6)=0.5 

XL(7)=0.6 

XL(8)=0.7 

XL(9)=1.0 

XL(10)=2.0 

XL(11)=3.0 

XL(12)=4.0 

XL(13)=7.7143 

XL(14)=10.0 

XL(15)=15.1429 

XL(16)=16.0 

XL(17)=17.0 

XL(18)=18.0 

XL(19)=19.0 

XL(20)=20.0 

XL(21)=20.1 
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XL(22)*20.2 
XL(23)=20.3 
XL(24)*20.4 
XL(25)=20.4167 
DO  102  N=l,25 
XL(N)  -  (L/20.)*XL(N)  -  L/2. 

102  CONTINUE 
BR(1)=0.0 
BR(2)=0.485 
BR(3)=0.658 
BR(4)=0.778 
BR(5)=0.871 
BR(6)=0.945 
BR(7)=1.010 
BR(8)=  1.060 
BR(9)=1.18 
BR(10)=1.41 
BR(11>=1 .57 
BR(12)=1.66 
BR(13)=1.67 
BR(14)=1.67 
BR(15)=1.67 
BR(16)=1.63 
BR(17)-1.37 
BR(18)=0.919 
BR(19)=0.448 
BR(20)=0. 195 
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BR(21)=0.188 
BR(22}=0.168 
BR(23)-0.132 
BR(24)=0.053 
BR(25)=0.0 
DO  104  K=l,25 
VEC0(K)=BR(K) 

VEC 1  (K)=XL(K)*BR(K) 
VEC2(K)=XL(K)*XL(K)*BR(K) 

104  CONTINUE 

CALL  TRAP(25,VECO,XL,EO) 

CALL  TRAP(25,  VEC  1  ,XL,E  1 ) 

CALL  TRAP(25,VEC2,XL,E2) 

C 

ALPHA=0.0 

ZD”ZDS+ALPHA*ZDB 

MD=MDS+ALPHA*MDB 

C  CALCULATING  THE  SS  PITCH  ANGLE  0O 
C  WITH  NEWTON  RAPHSON  METHOD 
Pl=  ZW*MD  -  MW*ZD 
P2=  XG*BUO*ZD 
P3=  ZG*BUO*ZD 
P4=  CD*(MD*E0  -  ZD*E1) 

WRITE(*,*)  P1,P2,P3,P4 
EPSI=. 00000001 
THETAO=0 
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IF(XGB.GT.O)  P4*-l  *P4 

18  DO  19  1=1,2000 
FT=FUNC(THET  AO) 

DFT=DFUNC(THETAO) 

DELT=FT/DFT 

THET  AO“THET  AO-DELT 
IF(ABS(DELT)-EPSI)  20,20,19 

19  CONTINUE 
FT=FUNC(THET  AO) 

20  WRITE(  V)  THETAO ,  FT 
C 

D  V =(MASS-ZWDOT)*(  I Y -MQDOTMMAS  S  ’XG+ZQDOT)  *  (MAS  S  *  XG+MWDOT) 
A11DV=(IY-MQDOT)*(ZW-2*CD*EO*U*TAN(THETAO)) 

&  +(MASS*XG+ZQDOT)*(MW+2*CD*E  1  *U*TAN(THETAO)) 

A!  2DV=(IY-MQDOT)*(MASS+ZQ+2*CD*El  *U*TAN(THETAO))+ 

&  (MASS*XG+ZQDOT) 

&  *(MQ-MASS*XG-MASS*ZG*U*TAN(THETA0)-2*CD*E2*U*TAN(THETA0» 
A 1 3  D  V=-(MAS  S  *  XG+ZQDOT)  *  WEIGHT 
BIDV  =(IY-MQDOT)*ZD+(MASS*XG+ZQDOT)*MD 
A21DV=(MASS-ZWDOT)*(MW+2*CD*El  *U*TAN(THETA0» 

&  +(MASS*XG+MWDOT)*(ZW-2*CD*EO*U*TAN(THETAO)) 

A22DV=(MASS-ZWDOT)*(MQ-MASS*XG-MASS*ZG*U*TAN(THETAO)-2*CD*E2*U 
&  *TAN(THETAO))+(MASS*XG+MWDOT)* 

&  (MASS+ZQ+2*CD*E  1  *U*TAN(THET  A0» 

A23  DV=-(MAS  S-ZWDOT)  *  WEIGHT 


B2DV  ^MASS-ZWIX)T)*MD+(MASS*XGfMWDOT)*ZD 
C 

A11=A11DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 
B1  =B1DV/DV 
B2  “B2DV/DV 
C 

EPS  =l.D-5 
ILMAX=1500 
C 

DO  1  1=1, ITC 
WRITE  (*,2001)  I,ITC 

TC=TCMIN+(I- 1  )*(TCMAX-TCMIN)/(ITC- 1 ) 

POLE=1.0/TC 
ALPHA3=4  0*POLE 
ALPHA2=6.0*POLE*  *2 
ALPHA1=4.0*POLE**3 
ALPHA0=  POLE**4 

K4=ALPHA0/((B  1  *  A2 1 -B2*  A1 1)*U0*  *4+(B  1  *  A23-B2*A13)*ZGB*U0**2) 

A2M=B1*U0**2 

A3M=B2*U0**2 

A0M=-(A1 1  +A22)*U0-ALPHA3 

B1M=B2*U0**2 
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B2M«(B2*A12-B  1  *A22)*U0**3 
B3M=(B  1  *A21-B2*A1 1)*U0**3 

B0M=(A1 1  *  A22-A2 1  *A12)*U0**2-A23  *ZGB-ALPHA2-B  1  *U0*U0*K4 

C1M=(B2*A11-B1*A21)*U0**3 

C2M=(B1*A23-B2*A13)*ZGB*U0**2 

C0M=(A13*A21-A23*A1 1)*ZGB*U(H-ALPHA1-(B2+B1  *A22-B2*A12)*K4*U0**3 
K2=C1M*B0M*A3M-B  1  M*C0M*A3M-C  1M*B3M*  AOM 
K2=K2/(C1M*B2M*A3M-B1M*C2M*A3M-C1M*B3M*A2M) 

K 1  =(C0M-C2M*K2)/C  1 M 
K3=(AOM“A2M*K2)/A3M 

DO  2  J=1,IU 

U-  UMIN+(  J- 1  )*(UMAX-UMIN)/(IU- 1 ) 

A(1.1H).0D0 

A(1,2>=0.0D0 

A(1,3)=1.0D0 

A(l,4)=0.0D0 

A(2, 1  )=-A13*(XGB*SIN(THETAO)-ZGB*COS(THETAO))+B  1  *U*U*K1 
A(2,2)=A1 1  *U  +B1*U*U*K2 
A(2,3)=A12*U  +B1*U*U*K3 
A(2,4)=  B1*U*U*K4 

A(3 , 1  )=- A23  *(XGB  *  S  IN(THET  A0)-ZGB  *  COS(THET  A0))+B2*  U*  U*K  I 

A(3,2)=  A21*U  +B2*U*U*K2 

A(3,3)=  A22*U  +B2*U*U*K3 

A(3,4)=  B2*U*U*K4 

A(4,l)=-  U 

A(4,2)s=  1.0D0 
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A(4,3)=  O.ODO 
A(4,4)=  O.ODO 
C 

CALL  RG(4,4,  A,WR,  WI,0,ZZZ,IV  1  ,FV  1  ,IERR) 
CALL  DSTABL(DEOS,WR,WLFREQ) 

C 

IF  (J.GT.l)  GO  TO  10 
DEOSOO  DEOS 
UOO  -  U 
LL=  0 
GO  TO  2 

10  DEOSNN  =  DEOS 
UNN  =  U 

PR=  DEOSNN*DEOSOO 
IF  (PR.GT.0.D0)  GO  TO  3 
LL  =  LL+1 

IF  (LL.GT.3)  STOP  1000 
IL  =  0 
UO  =  UOO 
UN  =  UNN 
DEOSO  =  DEOSOO 
DEOSN  =  DEOSNN 
6  UL  *  UO 
UR  -  UN 
DEOSL  =  DEOSO 
DEOSR  =  DEOSN 
C  U  =  (UL+UR)/2.D0 
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ALPHA  -  (DEOSL-DEOSR)/(UL-UR) 

U  =-  (DEOSL-ALPHA*UL)/ALPHA 

A(l,l>  =  0.0D0 

A(l,2)  =  O.ODO 

A(l,3)  =  1.0D0 

A(l,4)  -  O.ODO 

A(2,l)  =  A13*(XGB*SIN(THETAO)-ZGB*COS(THETAO))+Bl*U*U*Kl 
A(2,2)  =  A11HJ  +B1*U*U*K2 
A(2,3  )  =  A12*U  +B1*U*U*K3 
A(2,4  )=  B1*U*U*K4 

A(3,l)  =  A23*(XGB*SIN(THETAO)-ZGB*COS(THETAO))+B2*U*U*Kl 

A(3,2)  =  A21*U  +B2*U*U*K2 

A(3,3)  =  A22*U  +B2*U*U*K3 

A(3,4  )  =  B2*U*U*K4 

A(4,l)=-U 

A(4,2)  =  l.ODO 

A(4,3)  =  O.ODO 

A(4,4  )=  O.ODO 

CALL  RG(4,4AWR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

DEOSM  =  DEOS 
UM  =  U 

PRL  =  DEOSL*DEOSM 
PRR  =  DEOSR'DEOSM 
IF  (PRL.GT.O.DO)  GO  TO  5 
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■ 


J 


UO  *  UL 
UN  =  UM 
DEOSO  =  DEOSL 
DEOSN  =  DEOSM 
IL  =  DL+1 

IF  (EL.GT.ILMAX)  STOP  3100 
DIF  =  DABS^UL-UM) 

IF  (DIF.GT.EPS)  GO  TO  6 
U  =  UM 
GO  TO  4 

5  IF  (PRR.GT.0.D0)  STOP  3200 
UO  =  UM 
UN  =  UR 
DEOSO  =  DEOSM 
DEOSN  =  DEOSR 
IL  =  IL+1 

IF  (IL.GT.ILMAX)  STOP  3 100 
DIF  =  DABS(UM-UR) 

IF  (DIF.GT.EPS)  GO  TO  6 
U  =  UM 

4  LLL  =  10+LL 

WRITE  (LLL.*)  U/U0,TC*U0/L 
3  UOO  =  UNN 

DEOSOO  =  DEOSNN 
2  CONTINUE 
1  CONTINUE 
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2001  FORMAT  (215) 
END 


SUBROUTINE  DSTABL(DEOS,WR,WI,OMEGA) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  WR(4),WI(4) 

DEOS=-1.0D+20 
DO  1 1  =  1,4 

IF  (WR(I).LT.DEOS)  GO  TO  1 
DEOS  =  WR(I) 

U  =  I 

1  CONTINUE 
OMEGA  -  WI(IJ) 

OMEGA  =  DABS(OMEGA) 

RETURN 

END 

C 

SUBROUTINE  TRAP(N,A,B,OUT) 

C  NUMERICAL  INTEGRATION  ROUTINE  USING  THE  TRAPEZOIDAL  RULE 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DIMENSION  A(1),B(1) 

N1  =  N-l 
OUT  =  0.0 
DO  1 1  =  1,N1 

OUT1  =  0 . 5  *(  A(I)+ A(I+ 1 ))  *(B(I+ 1  )-B(T)) 

OUT  =  OUT+OUT1 
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1  CONTINUE 
RETURN 
END 

C  FUNCTIONS  USED  IN  NEWTON  RAPHSON  ROUTINE 
FUNCTION  FUNC(THETAO) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  P1,P2,P3,P4 

FUNC  =  P1*DTAN(THETA0)  +  P2*DCOS(THETAO)  +  P3*DSIN(THETA0) 

&  +  P4*(DTAN(THETA0))**2 

RETURN 

END 

FUNCTION  DFUNC(THETAO) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

COMMON  P1,P2,P3,P4 

DFUNC=P1*(1.DO/DCOS(THETAO))**2  -  P2*DSIN(THETA0)  + 

&  P3*DCOS(THETAO)  +  P4*2.DO*DTAN(THETAO)*(1.DO/DCOS(THETAO))**2.DO 
RETURN 
END 
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